function y = aft_GetVolume(hfile, volnum)
% AFT::GetVolume  - get one volume from a (multi-volume) dataset
%
% FORMAT:       y = obj.GetVolume([volnum])
%
% Input fields:
%
%       volnum      1x1 double index (default: 1)
%
% Output fields:
%
%       y           uninterpolated volume data (in source datatype)
%
% Note: this methods works with
%
%       amr       get AMR data (volnum discarded, can only be 1)
%       cmp       get component map with number volnum
%       ddt       get the corresponding volume (1 .. 12)
%       dmr       get according volume of DWI file
%       fmr       get according volume of STC file(s)
%       glm       get beta map with index volnum (scheme PxS)
%       hdr       get sub-volume of 4D image file (1st volume otherwise)
%       head      get sub-brick data from HEAD/BRIK object
%       map       get FMR based map (volnum discarded, can only be 1)
%       msk       get MSK data (volnum discarded, can only be 1)
%       vdw       get tensor imaging volume with number volnum
%       vmp       get according statistical map with number volnum
%       vmr       get VMR data (volnum discarded, can only be 1)
%       vtc       get functional volume with number volnum
%
% TYPES: AMR, AVA, CMP, DDT, DMR, FMR, GLM, HEAD, MAP, MSK, VDW, VMP, VMR, VTC

% Version:  v0.8a
% Build:    9102123
% Date:     Oct-21 2009, 11:23 PM CEST
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 1 || ...
    numel(hfile) ~= 1 || ...
   ~isBVQXfile(hfile, true)
    error( ...
        'BVQXfile:BadArgument', ...
        'Invalid call to ''%s''.', ...
        mfilename ...
    );
end
if nargin < 2 || ...
   ~isa(volnum, 'double') || ...
    numel(volnum) ~= 1 || ...
    isinf(volnum) || ...
    isnan(volnum) || ...
    volnum < 1
    volnum = 1;
else
    volnum = floor(volnum);
end

% get super-struct
ssc = bvqxfile_getscont(hfile.L);
ft = lower(ssc.S.Extensions{1});
if ~any(strcmp(ft, ...
    {'amr', 'ava', 'cmp', 'ddt', 'dmr', 'fmr', 'glm', 'hdr', ...
     'head', 'map', 'msk', 'vdw', 'vmp', 'vmr', 'vtc'}))
    error( ...
        'BVQXfile:BadArgument', ...
        'GetVolume not supported for this object type.' ...
    );
end

% depends on filetype
bc = ssc.C;
switch (ft)
    
    % AMR
    case {'amr'}
        
        % initialize data
        y = uint8([]);
        y(1:size(bc.Slice(1).AMRData, 1), 1:size(bc.Slice(1).AMRData, 2), 1:numel(bc.Slice)) = 0;
        
        % iterate over slices
        for sc = 1:numel(bc.Slice)
            y(:, :, sc) = bc.Slice(sc).AMRData(:, :);
        end
        
    % AVA
    case {'ava'}
        am = fieldnames(bc.Maps);
        while volnum > 1 && ...
           ~isempty(am)
            avm = bc.Maps.(am{1});
            if size(avm, 4) < volnum
                volnum = volnum - size(avm, 4);
                am(1) = [];
            else
                y = avm(:, :, :, volnum);
                return;
            end
        end
        if isempty(am)
            error( ...
                'BVQXfile:BadArgument', ...
                'Volume number out of bounds.' ...
            );
        end
        avm = bc.Maps.(am{1});
        y = avm(:, :, :, volnum);
        
    % CMP
    case {'cmp'}
        if volnum > numel(bc.Map)
            error( ...
                'BVQXfile:BadArgument', ...
                'Volume number out of bounds.' ...
            );
        end
        y = bc.Map(volnum).CMPData(:, :, :);
        
    % DDT
    case {'ddt'}
        if volnum > 12
            error( ...
                'BVQXfile:BadArgument', ...
                'Volume number out of bounds.' ...
            );
        end
        if ~istransio(bc.TensorEigenVs)
            y = squeeze(bc.TensorEigenVs(volnum, :, :, :));
        else
            y = resolve(bc.TensorEigenVs);
            y = squeeze(y(volnum, :, :, :));
        end
        
    % DMR
    case {'dmr'}
        if volnum > bc.NrOfVolumes
            error( ...
                'BVQXfile:BadArgument', ...
                'Volume number out of bounds.' ...
            );
        end
        
        % depends on storage format
        switch (bc.DataStorageFormat)
            case {2}
                y = squeeze(bc.DWIData(:, :, volnum, :));
            case {3}
                y = bc.DWIData(:, :, :, volnum);
            case {4}
                y = squeeze(bc.DWIData(volnum, :, :, :));
            otherwise
                error( ...
                    'BVQXfile:BadObject', ...
                    'DataStorageFormat unknown.' ...
                );
        end
        
    % FMR
    case {'fmr'}
        if volnum > bc.NrOfVolumes
            error( ...
                'BVQXfile:BadArgument', ...
                'Volume number out of bounds.' ...
            );
        end
        
        % data must be loaded
        if isempty(bc.Slice)
            error( ...
                'BVQXfile:BadObject', ...
                'STC data not loaded.' ...
            );
        end
        
        % depends on FileVersion and DataStorageFormat
        if bc.FileVersion < 5 || ...
            bc.DataStorageFormat == 1
            y = uint16([]);
            y(1:bc.ResolutionX, 1:bc.ResolutionY, 1:bc.NrOfSlices) = uint16(0);
            
            % iterate over slices
            for sc = 1:bc.NrOfSlices
                y(:, :, sc) = bc.Slice(sc).STCData(:, :, volnum);
            end
        else
            switch (bc.DataStorageFormat)
                case {2}
                    y = squeeze(bc.Slice.STCData(:, :, volnum, :));
                case {3}
                    y = bc.Slice.STCData(:, :, :, volnum);
                case {4}
                    y = squeeze(bc.Slice.STCData(volnum, :, :, :));
                otherwise
                    error( ...
                        'BVQXfile:BadObject', ...
                        'DataStorageFormat unknown.' ...
                    );
            end
        end
        
    % GLM
    case {'glm'}
        
        % what type of GLM
        if ~bc.ProjectTypeRFX
            nd = numel(size(bc.GLMData.BetaMaps));
            if volnum > size(bc.GLMData.BetaMaps, nd)
                error( ...
                    'BVQXfile:BadArgument', ...
                    'Volume number out of bounds.' ...
                );
            end
            if nd == 2
                y = bc.GLMData.BetaMaps(:, volnum);
            elseif nd == 4
                y = bc.GLMData.BetaMaps(:, :, :, volnum);
            else
                error( ...
                    'BVQXfile:BadObject', ...
                    'Unsupported GLMData size.' ...
                );
            end
        else
            nd = numel(size(bc.GLMData.Subject(1).BetaMaps));
            ns = size(bc.GLMData.Subject);
            nsp = size(bc.GLMData.Subject(1).BetaMaps, nd);
            if volnum > (ns * nsp)
                error( ...
                    'BVQXfile:BadArgument', ...
                    'Volume number out of bounds.' ...
                );
            end
            snum = floor((volnum - 1) / nsp) + 1;
            pnum = volnum - (snum - 1) * nsp;
            if nd == 2
                y = bc.GLMData.Subject(snum).BetaMaps(:, pnum);
            elseif nd == 4
                y = bc.GLMData.Subject(snum).BetaMaps(:, :, :, pnum);
            else
                error( ...
                    'BVQXfile:BadObject', ...
                    'Unsupported GLMData size.' ...
                );
            end
        end
        
    % HDR
    case {'hdr'}
        if volnum > size(bc.VoxelData, 4)
            error( ...
                'BVQXfile:BadArgument', ...
                'Volume number out of bounds.' ...
            );
        end
        y = bc.VoxelData(:, :, :, volnum);
        
    % HEAD
    case {'head'}
        if volnum > numel(bc.Brick)
            error( ...
                'BVQXfile:BadArgument', ...
                'Volume number out of bounds.' ...
            );
        end
        y = bc.Brick(volnum).Data(:, :, :);
        
    % MAP
    case {'map'}
        
        % initialize data
        y = single([]);
        y(1:bc.DimX, 1:bc.DimY, 1:numel(bc.Map)) = single(0);
        
        % iterate over slices
        for sc = 1:size(y, 3)
            y(:, :, sc) = bc.Map(sc).Data(:, :);
        end
        
    % MSK
    case {'msk'}
        y = bc.Mask(:, :, :);
        
    % VDW
    case {'vdw'}
        if volnum > bc.NrOfVolumes
            error( ...
                'BVQXfile:BadArgument', ...
                'Volume number out of bounds.' ...
            );
        end
        
        % make sure to load transio object
        t = bc.VDWData;
        if istransio(t)
            t = t(:, :, :, :);
        end
        y = squeeze(t(volnum, :, :, :));
        
    % VMP
    case {'vmp'}
        if volnum > numel(bc.Map)
            error( ...
                'BVQXfile:BadArgument', ...
                'Volume number out of bounds.' ...
            );
        end
        y = bc.Map(volnum).VMPData;
        
    % VMR
    case {'vmr'}
        if volnum == 1
            y = bc.VMRData;
        elseif ~isempty(bc.VMRData16)
            y = bc.VMRData16;
        else
            error( ...
                'BVQXfile:BadArgument', ...
                'Higher volume only valid if V16 loaded.' ...
            );
        end
        
    % VTC
    case {'vtc'}
        if volnum > bc.NrOfVolumes
            error( ...
                'BVQXfile:BadArgument', ...
                'Volume number out of bounds.' ...
            );
        end
        
        % make sure to load transio object
        t = bc.VTCData;
        if istransio(t)
            t = t(:, :, :, :);
        end
        y = squeeze(t(volnum, :, :, :));
end

% final check
if istransio(y)
    y = y(:, :, :);
end
